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Abstract 

We present a phenomenological model for the initial conditions needed in a one- 
fluid hydrodynamical description of ultrarelativistic nuclear collisions at CERN- 
SPS. The basic ingredient is the parametrization of the baryon stopping, i.e. the 
rapidity distribution, as a function of the thickness of the nuclei. We apply the 
model to S + S and Pb + Pb collisions and find after hydrodynamical evolution 
reasonable agreement with the data. 
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1 Introduction 


One of the main motivations of heavy ion experiments is the investigation of the 
nuclear equation of state away from the ground state. However, it is very difficult 
to deduce from the final multi-particle state, i.e. from the experimental particle 
distributions, model independent statements about the equation of state. It is 
not even granted, that the time and volume available are large enough to create a 
form of matter which can be described using (infinite volume) thermodynamics. 
Nevertheless, the only possibility to obtain information on the equation of state 
at large temperatures and densities is via dynamical models of these collisions. 

One popular approach to the dynamics of heavy ion collisions is the use of one- 
fluid hydrodynamics ill- Hydrodynamical simulations at energies below 10 A 
GeV (BEVALAC and SIS) usually start with the approaching nuclei before their 
initial impact and include the initial compression. In such a treatment the nuclei 
fuse to a single fluid, implying, at zero impact parameter, a complete stopping 
of equal size nuclei. At higher energies, as at the CERN-SPS, RHIC and LHC, 
this is not justified, and one must be able to incorporate nuclear transparency in 
the description. Instead of trying to describe the production and equilibration 
within hydrodynamics, one starts the calculation from initial conditions which 
specify the hydrodynamic state of the system at initial time t q. Initial conditions 
parametrize the production and equilibration dynamics. 

The aim of this article is to present a parametrization of the initial baryon 
stopping, i.e. initial baryon rapidity distribution, in terms of the thickness of the 
incoming nuclei. The parametrization should be valid for all symmetric collisions 
at a given cms energy. Here we concentrate on CERN-SPS collisions at yfs ~ 20 
GeV per nucleon. After having found such a parametrization there are only a few 
additional parameters left like, e.g. the initial time. This makes the determination 
of the initial state from experimental data, usually done by trial and error, less 
arbitrary. 

2 Baryon Stopping 

The experimental net baryon rapidity distribution dN B /dy is in general a function 
of y/s and the mass numbers A and B of the colliding nuclei. The mass numbers 
determine the mean number of interactions which an individual nucleon suffers. 
The rapidity distribution of final free baryons, dN B /dy {yfs, A, B), is measured 
in several experiments. However, since we want to study the evolution of the 
colliding system after its formation at time to, we need the knowledge of dN B /dy 
at to- This information can only be obtained by models, which describe the 
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Figure 1: a) Rapidity distribution of protons in p + p collisions for three different 
energies normalized to their maximum: Experimental data: ’x’ y/s = 4.9 GeV 0, ’o’ 
= 6.9 GeV 0 , ’+’ Vs = 27.5 GeV The line corresponds to the parametriza- 
tion (Jl|) described in the text. 

b) Parametrization of the rapidity distribution of the initial baryons for various num¬ 
bers of average nucleon-nucleon collisions (n) normalized to the maximum. Going 
from the edge to the center the lines correspond to (n) =0.1 (dotted), (n) = 1 (full), 
(n) = 2 (dashed-dotted), (n) = 4 (dashed) and (n) = 8 (dotted). The line for (n) = 1 
is the same as the parametrization in Fig. [jja. 


initial production stage of the nuclear collision well. Cascade type models may 
in principle provide this information. Our approach is much simpler. All physics 
of baryon stopping is included in the parametrization of dN B /dy with the final 
justification coming from the comparison to the experiment. 

To motivate our parametrization we start from what is known in p + p col¬ 
lisions. In order to compare various collisions at different energies, we use the 
scaled rapidity x y = y/y m ax, where y is measured in the cms frame of the col¬ 
liding protons. In Fig. jl]a we show the rapidity distribution dN/dy of protons 
for three values of y/s, plotted against x y and normalized to 1 at the maximum. 
Going from lower to higher energies, the central region (y ~ 0) becomes more 
depleted indicating increasing transparency of colliding nucleons with increasing 
y/s. The maximum, however, seems to stay fixed around x y = 0.85 independent 
of y/s. 

In Fig. [|a we also show a parametrization of the dN/dy as we interpolate 
the data to y/s = 20 GeV as appropriate for CERN-SPS nuclear collisions. The 
following consideration went into the parametrization. We compose the distribu¬ 
tion out of two parts, one dominant in the target and the other in the projectile 
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fragmentation region. The rise in dN/dy going from the center to forward (or 
backward) regions is roughly exponential. Therefore we start out with two expo¬ 
nential functions of third order polynomials: 


dN 

dy 


dN Pm J, dN Tar 

' \ x y) + 


dy 


dy 


(x y ) = 


Cp ro j exp(a.T 3 + bxy + cx y ) + Ct ar exp(-ax 3 + bx'l 


(1) 



(1 



The factor (1 — x y ) is chosen to cut off the distribution at the phase space bound¬ 
ary. The fit in Fig. |ja is given by a = 1.5, b = 0, c = 3.0, and Cp ro j = Cx ar with 
values of C’s determined via the normalization. 

We now assume that the initial baryon rapidity distribution dN-g/dy in a 
nuclear collision has the same functional form as Eq. but the parameters a, b 
and c depend on the nuclear thickness a participating beam (target) nucleon sees 
on its way through the target (beam) nucleus. The task is to find this dependence. 
The thickness function is defined as 

Ta(p) = J dz p A {z,p) , (2) 

where pa(t) is the nuclear density for a nucleus of mass number A, z is the 
longitudinal and p the transverse variable: f = ( z,p ). We take a Woods-Saxon 
parametrization for the nuclear density 

PAif) = - jTpj §-T7 f— T ) (3) 

exp[(|F| - R A )/a R \ + 1 

with 

Ra = 1.12 fm x A 1/3 - 0.86 frn x A" 1/3 (4) 

a R = 0.54 fm, and po = 0.17 fm -3 |7|. 

We expect that the stopping depends not only on the thickness but also on an 
interaction strength, like the inelastic cross section. Then the natural quantity 
for the parametrization of stopping is 


(n){p) = a pp T A (p) , (5) 

where we use a pp = 32 mb, the total inelastic cross section for p + p collisions at 
SPS energy. In the Glauber model || (n) (p) is the average number of interactions 
suffered by a nucleon colliding with a nucleus at impact parameter b = p. The 
value of (j pp is not important since (n) will be multiplied by adjustable parame¬ 
ters. We should like to emphasize that ( n ) should be considered as a reasonable 
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quantity to describe the strength of stopping of a nucleon and it need not be 
interpreted as the number of interactions except at the limit of small (n). In 
particular, for (n) = 1 we should have the same rapidity distribution as in p + p 
collisions. 

We impose the following limits. Since \x y \ < 1 the parameter a becomes 
important at the edge of the phase space. For Ta —► oo we expect a to vanish 
because the nucleons are shifted gradually to x y = 0. The limit a — > 0 is achieved 
by a oc 1 /Ta- More generally we could assume power behaviour but the inverse 
dependence turns out to be sufficient. 

In Eq. <® the actual stopping is expressed through the parameter b which is 
zero for (n) = 1 and negative for (n) > 1, and should dominate for growing (n). 
We use a simple ansatz b = /3 S (1 — (n)) = (3 S (1 — <j pp Ta)- The prefactor f} s is a 
fit parameter characterizing the stopping. From fits to the experimental particle 
spectra we obtain (3 S = 2.25. One has to keep in mind that the deduced value of 
(3 S depends on the time evolution from initial to final distributions and therefore 
on the equation of state which is used. Finally, the parameter c could also be 
chosen to depend on the thickness. However, it turns out that we can obtain the 
needed rapidity shift as a function of nuclear thickness from the ( n ) dependence 
of b((n)) alone and leave c as a constant. 

We summarize the functional dependence of the parameters a, b and c on Ta- 

a(T A ) = 1.5 (cjppTA) -1 

b(T A ) = A, (1 - a pp T A ) 

c(T a ) = 3.0. (6) 

These functions @ together with © form our phenomenological ansatz with the 
values of the parameters determined mainly from the S + S data. In Fig. 0b we 
show, with /3 S = 2.25 the resulting dN-g/dy for several values of (n) normalized 
as in Fig la. For (n) = 1 the curve is the same as in Fig. [Ija. The dependence 
of dNft/dy on (n) shows the following limits. For (n) —> 0 the parameter a((n)) 
dominates and leads to a sharp rise at \x y \ < 1 which is cut down by the phase 
space factor at the boundaries. Therefore dN-g/dy turns to a double 5-function 
peaked at target and projectile rapidity. For (n) —> oo, the parameter b((n)) 
dominates and results in a Gaussian function for dN^/dy which narrows with 
increasing (n). 

Since we obtain the distribution as a function of transverse variable dN^(p)/dy 
it will be normalized to the number of nucleons per unit transverse area at p. To 
obtain the final local density we must know the dependence between rapidity and 
longitudinal distance z at the initial time to- We do this by specifying the initial 


5 


velocity profile by choosing a linear relation between rapidity and z: 

y{p, z) = k{p)z , v z (p, z) = tanh [k{p)z] . (7) 

The local density ps{z,p) is obtained from the dN^ip)/dy by multiplying with 
( dy/dz )/coshy where coshy is the Lorentz contraction factor between the local 
rest frame and the fixed cms frame. We will specify the initial flow velocity, in 
particular the p dependence of n(p), after discussing the initial energy distribu¬ 
tion. 


3 The Initial Energy Flow and Velocity Pro¬ 
file 


The numerical solution of the hydro dynamical equations is determined on a 2 + 1 
dimensional grid in the cylindrical coordinate system r = (z,p) [[|. The frame 
which we use is always the cms frame of the participating nucleons. 

Since we assume azimuthal symmetry the simulation is strictly valid only for 
zero impact parameter collisions. However, even the triggering to 5% of most 
central events in the experiments corresponds to a considerable range in impact 
parameter. It turns out that if we incorporate the full energy available in the 
impact parameter zero collisions we slightly overshoot the experimental spectra. 
To account for the experimental impact parameter averaging, we use effective 
nuclear sizes, i.e. we replace the incoming nuclei A. B by A eS = S,A and B c{i = £B 
and fix the geometry in terms of effective mass numbers and B cS . The 
value of £ is close to one. Our numerical algorithm for the calculation of final 
particle spectra leads to a few per cent loss of baryon number and energy ||]. We 
compensate also for these losses when fixing the value of £ to be 0.9 in symmetric 
collisions. 

For the hydrodynamical calculation we need also the local energy density e(r) 
and the initial flow velocity at the initial time to- For the energy density we follow 
a similar approach as for the baryon density. We first parametrize the energy per 
unit transverse area and unit rapidity using a Gaussian distribution in rapidity 


dE(p) 

dy 


C £ exp 


~(y - yo f 

2 (J I 


1 iy/ymax) 


( 8 ) 


where the width <r e and the normalization C £ depend on the transverse coordinate 
p. The last factor is the same phase space cut-off as imposed for the baryons 
©■ Eq. (|8|) contains only the non-baryonic energy to which the contribution 
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connected with the net baryon density will be added. The parametrization has 
three unknowns (at each p), but only one, e.g. the width a £ may be freely chosen. 
The normalization C e is given by energy conservation, imposed per unit transverse 
area at each value of p. The value yo is the center of mass rapidity for nucleons 
at p. For symmetric collision yo = 0 in the cms frame. 

It turns out that the experimental pion rapidity distribution is well reproduced 
when the value of a £ correlates with the baryon stopping. A larger stopping of 
baryon number results in a larger stopping of energy, too. Therefore the width 
of the energy distribution in rapidity space decreases with increasing nuclear 
stopping. We describe this with an ansatz 

&s(p) 7 tt, 7 / \i a e 5 (9) 

[Cpp Ta \p) opp Tb ( P ) ] 

where a pp is included in order to make the denominator dimensionless. The 
numbers appearing in (|j) are determined from the fits to the symmetric collision 
systems: c £ = 0.75, a £ = 0.13. 

Next we add to the energy density Eq. ®) the energy density connected with 
the net baryon number. This additional energy is important at the edge of phase 
space and ensures that we take into account the rest mass and the thermal energy 
of every leading nucleon. We write the coordinate space density in the form 


e(y, z) 


dE(p) dy(p, z) 
dy dz 


+ (Pb(p,z), 


( 10 ) 


where £ is the energy carried by the baryon. We take £ = 1.3 GeV ~ m nuc + E t h 
where m nuc is the rest mass of a nucleon and E t h its thermal energy. Since pb is 
the net baryon distribution, the energy of baryon-antibaryon pairs is included in 
the first term. 

The expansion of the matter is very sensitive on the initial velocity profile. 
Since we consider zero impact parameter collisions with cylindrical symmetry 
only, we do not expect significant collective motion in transverse direction initially, 
and take the velocity in the transverse direction at to to vanish, i.e., v p (to,f) = 0. 

In the case of the Bjorken model [ 10 1 , the scaling ansatz for longitudinal 
velocity is v z = z/t and initial conditions are usually defined on fixed proper time 
To- At finite, albeit high collision energies, the longitudinal extent of the system 
is finite and the scaling assumption must break in the fragmentation regions. 
Since we perform the numerical calculations in the center-of-mass frame of the 
participating nucleons we specify the initial condition at a fixed time to dr this 
frame. On this line the scaling velocity is z/t$. Since scaling cannot hold, we have, 
for convenience, chosen the linear ^-dependence for rapidity instead of velocity 
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Figure 2: Local baryon density pB (upper) and ype (lower) in the z-p plane as 
they result from the parametrization. The initial distribution is for S + S collisions 
(left) and Pb + Pb (right). 


as given in Eq. ©• For small 2 this ansatz can approach the velocity profile of 
the scaling solution. The proportionality constant is now k instead of 1/to of 
the Bjorken model. We define Tq® = 1 / k(p = 0) as a parameter which can be 
regarded as an equilibration time scale in the same way as tq in the scaling case. 

We determine the p dependence of n{p) through the longitudinal extension 
of the system, zl{p), at different values of p. The p dependence arises from 
the variation of the nuclear thickness in the transverse plane. The difference of 
the longitudinal extension between the center and the edge is of the order of 
A 2 = zl(p = 0 ) — zl{p = Ra) ~ Ra/i, where 7 = cosh( jm) is the Lorentz 
contraction of the nucleus in the fireball rest frame. This form follows from 
nuclear geometry and the assumption that formation and equilibration times are 
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independent of p. We define 


zl(p) 


"0 + 


R\~p 2 


7 


( 11 ) 


where zq is the longitudinal extension of the matter at equilibration time in the 
thin disk limit. Instead of zq we use as an independent parameter with the 
assumption that the forward (backward) edge of the initial fireball, zl(p), moves 
with the target (projectile) rapidity y cms , which, from Eq. i) at p = 0 gives 

Zo = r 0 Perns ■ ( 12 ) 

7 

Combining the above definitions one gets 


K(p) 


Veins 

zl{p) ' 


(13) 


Even though the initial distributions approach zero in transverse direction as 
a consequence of the finite size of the nuclei (cf. Eq. (§|)), it is still practical for the 
numerical work to make them zero beyond some value of p. As a criterium we use 
the mean value of interactions, (n) and set the initial distributions to zero when 
(n) < 0.5. It is clear that even so the matter at the edge is not dense enough for 
thermalization but for a finite system we are forced to extend the calculation at 
the surface to densities were the hydrodynamics cannot be justified. 

We show in Fig. ||| the resulting initial baryon density distributions for S + S 
and Pb + Pb. The upper frames give the densities in the local comoving frame 
and the lower the spatial density in the overall rest frame of the collision. The 
difference comes from the 7 factor of the Lorentz contracted fluid cell. The figure 
clearly shows the increase in transparency with radius. We think this is more 
reasonable than using the same baryon distribution at the edge as at the center. 
In this figure the cut in (n) leads to a sharp discontinuity in radial directions 
near the maximum. However, this does not influence the dynamics, because the 
pressure gradients are given by the local densities, which are small and smooth 
in this region when going to larger p as seen in the two upper figures. 

In Tab. [l] we give the values of parameters which characterize the initial state. 
An interesting point is that Tq® which is a parameter adjusted to each collision 
system turns out to be rather similar in Pb + Pb and S + S collisions. This 
means that the equilibration time seems to be independent of the thickness of the 
nuclei. 
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Table 1: Summary of parameters characterizing the collision when a QGP equation 
of state is used. £ and Tq® are parameters adjusted to the data and explained in the 
text. The rest follows from the initial parametrization. tf is the central freeze-out 
time and (v p ) the mean radial velocity in the rapidity range of \y z < 0.25. 


collision 

s + s 

Pb + Pb 

£ 

0.9 

0.9 

Tg ff (fm/c) 

1.2 

1.3 

zl{0 ) (fm) 

3.5 

3.6 

cent, e (GeV/fm 3 ) 

7.1 

16.7 

cent, pb (fm -3 ) 

0.59 

1.94 

cent. T (MeV) 

244 

300 

R(£A) (fm) 

3.15 

6.25 

t f (fm/c) 

6.1 

11.4 

{ v p){\Uz\ < 0.25) 

0.28 

0.34 


4 Results 


We illustrate the use of the parametrization of initial conditions discussed above 
by comparing with experimental data from the NA35 [111 and NA49 || collab¬ 
orations. We solve the hydrodynamical evolution as described in ||, where also 
the treatment of the freeze-out is explained. We use isospin symmetry when cal¬ 
culating the hadron spectra. The isospin corrections would somewhat reduce the 
calculated p — p spectrum in the heavy systems (e.g. Pb + Pb) but our choice 
£ = 0.9 probably overestimates the amount of spectator nucleons. The freeze-out 
happens at a constant energy density of ey = 0.15 GeV/fm 3 leading to an average 
freeze-out temperature of Tf ~ 140 MeV. In [3j] we also explain the construction 
of the equation of state we use. Here we show results only for an equation of state 
with phase transition to QGP at T c = 165 MeV, labeled as EOS A in 

The stopping shows up mainly in the rapidity distributions. We show in Fig. |] 
the distributions for negatives and proton-antiproton difference in a central S + 
S collision. In the p — p rapidity distribution target and projectile fragmentation 
regions exhibit clear maxima, which are nicely reproduced by the calculation. The 
normalization is slightly adjusted by the £ parameter. The rapidity distribution 
of the negatives fits very well. Also the transverse momentum distributions of 
both negatives and p — p are reproduced reasonably well. 


10 










Figure 3: Comparison of experimental data from S + S collisions with final hadron 
spectra calculated using an equation of state (EOS A) with phase transition at T c = 
165 MeV. The data are taken by the NA35 collaboration |ll|. 


Fig. ^ shows the result for Pb + Pb collisions compared with preliminary 
data from the NA49 collaboration ||. There are small deviations in the proton 
distribution but the negative particle rapidity distribution is very well reproduced. 
Both collisions are calculated with nearly the same parametrization, the difference 
coming from the thickness function Tq (p) only. The only slight difference is the 
7 g ff = 1.2 for S + S and Tg ff = 1.3 for Pb + Pb. Thus the parametrization 
reproduces well the stopping as a function of nuclear size. The result for p — p 
show larger deviation with faster protons in the experimental spectrum than in the 
calculated result. At least part of this must come from the finite impact parameter 
range included in the measurement and leading to a larger average transparency. 
In the calculation we take this into account only in the normalization through 
the parameter £ with practically no change in the shape of the spectrum. 

The gross features of the transverse momentum data of the negative particles 
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Figure 4: Comparison of experimental data from Pb = Pb collisions with final 
hadron spectra calculated using an equation of state (EOS A) with phase transition 
at T c = 165 MeV. The data are taken by the NA49 collaboration ||] and should 
be considered as preliminary. The mr-spectra of negatives and p — p are measured 
in rapidity intervals of width 0.5 centered at 2.9 (p — p only), 3.4, 3.9, 4.4, 4.9, 
5.4. For clarity the data and the calculation are successively scaled down by 10 -n 
(n = 0,1,2,...). 


are again well reproduced except for the very forward (and backward) region 
were the calculated slope is steeper than the data. The particle density at this 
edge of the phase space, however, is becoming so small that the hydrodynamics 
with relatively strong longitudinal flow can lead to an artificially large transverse 
cooling. Thus we should not expect that hydrodynamics can describe the far 
edges of fragmentation regions well. 

This problem is more pronounced in the p — p transverse distributions which 
indicate that even in the central rapidity region the experimental transverse flow 
is somewhat stronger. A possible explanations is that the equation of state is 
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too soft. Later freeze-out time might also enhance the flow effect on protons as 
compared to pions. 

In Tab. [I] we also give numbers which characterize the freeze-out. The life¬ 
time tf of the fireball, measured at the center, scales as expected with the size of 
the fireball. We see a doubling of the lifetime going from S + S to Pb + Pb in 
accordance with the doubled radius of the lead nucleus. We further give the av¬ 
erage (radial) transverse velocity (v p ) of the fluid cells at freeze-out with rapidity 
\y z \ < 0.25. The result may be compared with an analysis of NA44 [ 12 1 . They 
extracted from a fit to hadron data a mean transverse velocity. The result is ]l2|] 
( v p ) = 0.24 ± 0.10 for S + S and ( v p ) = 0.36 ± 0.14 for Pb + Pb in agreement 
with our studies. However, the NA44 obtained a higher freeze-out temperature 
for Pb + Pb, Tf = 167 ± 13 MeV (for S + S Tf = 142 ± 5 MeV). This explains 
why NA44 gets agreement with ( v p ) = 0.36 while we slightly underestimate the 
slope in the central p — p transverse momentum spectra. 

Finally we compare our results with the related work of Schlei et al. [ 23 . (See 
also |l3]] in this Volume, where a different parameter set for initial conditions is 
used.) The main difference is in the geometry of the initial fireball. The equations 
of state differ slightly. Both exhibit a phase transition to QGP, ours with T c = 165 
MeV vs. 200 MeV in j|]. In addition, the baryonic pressure in [|[ is omitted. 

The basic differences may be summarized in the following way. The calculation 
in Q starts from a small longitudinal extension, indicating large compression, and 
resulting in larger initial energy density as compared to ours. At the center of the 
fireball for S+S collision e* = 13.0 GeV/fm 3 in (2j vs. 7.1 GeV/fm 3 in our case. 
For Pb+Pb the numbers are 20.4 GeV/fm 3 and 16.7 GeV/fm 3 , respectively. At 
first glance this looks like a difference in stopping but a closer inspection shows 
that it is related more closely to the initial volume than to the initial velocity 
distribution. This can be interpreted as a difference in initial time, our calculation 
corresponding to a later initial time. If the energy stopping is characterized with 
the fraction of thermal energy, the numbers are very similar: 0.43 (0.64) in |2) vs. 
our 0.45 (0.66) for S+S (Pb+Pb). The difference in initial volume is related to the 
longitudinal extent and our larger initial volume can be interpreted as starting 
the calculation at later initial time. This shows up also as shorter lifetimes in our 
calculation, 6.1 fm/c vs. 6.9 fm/c in [|| for S+S and 11.4fm/c vs. 13.5 fm/c in (H 
for Pb+Pb collisions. 
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5 Conclusions 


A model for initial conditions was developed with a parametrization of baryon 
stopping in terms of the nuclear thickness as the basic input. The same parametri¬ 
zation gives a very satisfactory description of the basic features of hadron spectra 
both in S + S and Pb + Pb collisions. We think that this is an improvement 
compared to earlier approaches !,! where the transverse dependence of stopping 
was not taken into account. 
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